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ABSTRACT 

Nearby gamma-ray bursts (GRBs) are likely to have represented a significant threat to life 
on the Earth. Recent observations suggest that a significant source of such bursts is compact 
binary mergers in globular clusters. This link between globular clusters and GRBs offers the 
possibility to find time intervals in the past with higher probabilities of a nearby burst, by 
tracing globular cluster orbits back in time. Here we show that the expected flux from such 
bursts is not flat over the past 550 Myr but rather exhibits three broad peaks, at 70, 180 and 
340 Myr ago. The main source for nearby GRBs for all three time intervals is the globular 
cluster 47Tuc, a consequence of its large mass and high stellar encounter rate, as well as the 
fact that it is one of the globular clusters which comes quite close to the Sun. Mass extinction 
events indeed coincide with all three time intervals found in this study, although a chance 
coincidence is quite likely. Nevertheless, the identified time intervals can be used as a guide 
to search for specific signatures of GRBs in the geological record around these times. 

Key words: globular clusters: general - Gamma-ray burst: general - Galaxy: kinematics and 
dynamics - Astrobiology 



1 INTRODUCTION 

Globular clusters, densely packed groups of old stars, can effi- 
ciently produce close stellar binaries by dynamical interactions of 
their member stars. Examples for such dynamically formed bina- 
ries include low-mass X-ray binaries (e.g.|Clark|1975||Katz|1975| >, 
cataclysmic variables (Pooley & Hut 2006) and milli-second pul- 
sars (msPSRs, Ransom 2008 , Abdo et al. 2010 ). The most extreme 
binaries found in globular clusters consist of two neutron stars \An-\ 
|derson et al.|1990) . Mergers of such binaries are believed to be the 
central engine of short gamma-ray bursts (GRBs) (Grindl ay et al.| 
|2006||Dado et al.|2009||Lee et al.|20T0| >, that produce brief, intense 
flashes of ionising radiation. In contrast to short bursts, long bursts 
are believed to originate from the death of short-lived massive stars 
(see |Gehrels et al.|2009| for a review on long and short bursts). It 
has been argued that the rate of short GRBs in the local universe is 
dominated by the merger of neutron star binaries formed in globular 
clusters jSalvaterra et al.|2008||Guetta & Stella|2 009). A link be- 
tween globular clusters and short GRBs is further supported by the 
presence of a short GRB remnant candidate in the Galactic globu- 
lar cluster Terzan 5 ( Domainko 2011a), observed in the very-high 
energy gamma-ray ( Abramowski et al. 201 1 2013]), X-ray i Eger et 



al. 



2010 



al 



|Eger & Doma inko 2012 ) and radio wave band ( Clapson et 



2011 1. Additional evidence for the GRB - globular cluster con- 
nection comes from spatial offsets of short GRBs from their host 
galaxies ( |Berger|2010l |Salvaterra et al.|2010[|Church et al.|2011) 
and the redshift distribution of such events (Salvaterr a et al.|2008| 
IGuetta & Stella|2009| >. 



Since globular clusters follow well-defined orbits around the 
Galaxy (Domainko 201 lb), their coupling with GRBs allows us to 
examine the long-standing question of the past history of gamma- 
ray flux on the Earth. (A similar approach for supernovae explod- 
ing in star clusters has been used in Svensmark (2012)). Numer- 
ous studies have shown that gamma rays from supernovae (SNe) 
or GRBs could, in principle, have had a significant impact on the 
Earth's atmosphere and biosphere, potentially even contributing to 
mass extinctions (see Thorsett (1995); Scalo & Wheeler (2002); 
|Melott et al.l ( [2004l >; [Thomas et al.|p005l > for the affect of GRBs 
in general, and |Dar et aLl (l998l; Mel ott & Thomas| ( |201l] l for the 
affect of merger-induced bursts). However, demonstrating that SNe 
or GRBs may in fact have played some role first requires identify- 
ing that sources could have come near enough to the Earth at some 
point. Some previous studies have attempted to make a connection 
between the solar motion relative to the Galactic plane or spiral 
arms, on the assumption that the gamma ray flux incident on the 
Earth is larger in these regions of enhanced massive star formation 
rate and/or increased stellar density (see Bailer-Iones ( 2009]) for a 
review). However, a recent study shows that the flux from these 
sources as modulated by the plausible solar motion over the past 
550 Myr has a poor correlation with the variation of the extinction 
rate on the Earth (Feng & Bailer-Iones, submitted). 

Indeed, it seems that asu'onomical phenomena alone are un- 
likely to be the dominant driver of biological evolution or the cause 
of all (or even most) mass extinctions. Nonetheless, if a GRB were 
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to explode near to the Earth, its consequences could be catastrophic, 
and globular clusters are presumably a significant source of GRBs. 

The goal of this paper is to reconstruct the orbits of globular 
clusters relative to the Sun in order to calculate the GRB flux at 
the Earth as a function of time, and thereby to identify potential 
candidate clusters. The data for this orbital reconstruction comes 
from the positions, distance, proper motion and radial velocity cat- 
alogues of globular clusters of JDinescu et al.|1997|[T999a|bl[^03l >, 
from which we obtain the current Galactic coordinates and space 
velocities. By sampling over the (often significant) uncertainties in 
the reconstructed orbits of the globular clusters and the Sun, we 
infer the expected GRB flux as a function of time. This allows us 
to identify the most probable intervals in the Earth's history of a 
significantly increased gamma ray flux, which may (or may not) be 
associated with times of higher extinction rate. 

In section |2""fl we describe the orbital reconstruction method, 
and in section |2~2| we explain how we derive from this the proba- 
bility distribution over the past cluster-Sun separation and the ex- 
pected gamma ray flux at the Earth. This takes into account the 
different GRB rates in the clusters, which is derived in section |23l 
We give our results in section [3] where we also identify some past 
extinction events. We conclude in section|4]with an outlook on how 
to further this work. 



2 METHODS 

2.1 Reconstructing Galactic orbits 

We trace the orbits of the Sun and the globular clusters back in 
time by integrating the equations of motion through the Galactic 
potential. In a purely gravitational system there is no dissipation 
of energy, so the dynamics are reversible. We adopt an analytic, 
three component, axisymmetric potential, O, comprising the Galac- 
tic bulge, halo and disk 

<&(R,z) = <b b + ® h + ® d . (1) 
The bulge and halo are represented with a Plummer distribution 

~ GM ''J' 

= = (2) 



in which the characteristic length scales are b b = 0.35 kpc for the 
bulge and b h = 24.0 kpc for the halo, and the bulge and halo masses 
are M b = 1.40 x 10 10 M o and M h = 6.98 x 10" M G respectively. 
R is the radial coordinate perpendicular to the axis, and z is the 
distance from the Galactic plane. For the disk we use the potential 
from|Miyamoto & Nagai|(l975) 

-GM d 



4>,, 



R 2 + (a d + yjz 2 +bij 



(3) 



with the values M d 
a d = 3.55 kpc and b d = 0.25 kpc for the scale length and scale 
height of the disk, respectively (after Garci'a-Sanchez e t al.| ( |206"T) ). 
The integration is performed numerically from the present back to 
550 Myr BP (before present). This time limit is chosen because it 
corresponds to the beginning of the Phanerozoic eon, a time from 
which the fossil record becomes more indicative of biodiversity 
variations. The globular clusters (and Sun) are treated as massless. 

The initial conditions for the integration are the current phase 
space coordinates (three position and three velocity components) of 
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Figure 1. Samples of the orbit of 47 Tuc relative to the Sun to show how 
their separation varies over time. The variance arises from sampling the un- 
certainty in the current phase space coordinates of both the globular cluster 
and the Sun, and integrating each back in time through the Galactic poten- 
tial. 



the globular clusters (and Sun). These of course have significant un- 
certainties, each represented as a Gaussian with known mean (mea- 
sured coordinate) and standard deviation (estimated uncertainty). 
These come from Dana Casetti-Dinescu's catalogue for globular 
cluster's three-dimensional space velocities (20 12 version^] for the 
globular clusters, and from Hipparcos data by (Dehnen & Binney 
1998 ) for the Sun. We further use a distance of the Sun to the Galac- 
tic center obtained from astrometric and spectroscopic observations 
of the stars near the supermassive black hole of the Galaxy ( |Eisen-| 
|hauer et al.|2 003 ) and the displacement of the Sun from the Galactic 
plane is calculated from the photometric observations of classical 
Cepheids by Ma jaess et al.|p003| ). Rather than just performing a 
single integration for each object (cluster or Sun), we Monte Carlo 
sample its initial conditions from the uncertainty distribution in or- 
der to build up a large sample of orbits. Figure^shows an example 
of such sample orbits for one globular cluster, 47 Tuc, by plotting 
the distance of the cluster from the Sun over time. (We sample over 
the possible orbits of the Sun too.) We do not take into account 
the (possibly significant) uncertainties in the Galactic potential. In 
principle we could adopt an uncertainty model for these parameters 
and marginalize over them also. But we choose to omit this in this 
first investigation. 

Finally we have to note that compact binaries may be ejected 
from their parent cluster before they merge and produce a GRB 
(e.g. |Phinney & Sigurdsson|1991[|Ivanova et al.|20"08"} . This effect 
will smear out the distribution of compact binaries around the pro- 
ducing cluster. The typical escape velocities for massive globular 
clusters are about 50kms -1 , which is comparable to the present 
uncertainties of the globular cluster velocity. Although over time 
the orbit of the ejected binary could deviate considerably from its 
parent cluster, the uncertainty in its orbit is comparable to the uncer- 
tainty for its parent cluster, which we take into account. We there- 
fore choose to omit the issue of ejected GRB progenitors for this 
first investigation. Furthermore, more massive clusters are better 
able to retain their binaries, and these are the clusters that preferen- 
tially produce GRBs (see Sec.|2.3[l. 
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2.2 The probability distribution over globular cluster 
distances and the expected GRB flux at the Earth 

For a given globular cluster, c, we convert the set of (thousands 
of) relative orbits into a two-dimensional density distribution over 
time, t, and separation, r, using kernel density estimation. We in- 
terpret the resulting distribution as a probability distribution of the 
Sun-cluster separation over time, f c (r, f), which is normalized such 
that f r f c (r,i)dr = 1 for all t and for each cluster. This is shown 
in Figure [2] for 47Tuc, in which the probability density is plotted 
as a grey scale. At any given time, the darker the band, the more 
concentrated the probability is around a smaller range of distances. 
The width of the distribution at any time is determined by how the 
uncertainties in the present coordinates of both globular cluster and 
Sun propagate back in time. The density estimates for some other 
globular clusters are shown in Figures[5J|7] 

The flux of a gamma ray burst at the Sun is proportional to 
l/r 2 . Multiplying f c (r,t) by l/r 2 , and assuming that gamma ray 
bursts occur at random time^] we get a 2D distribution which is 
proportional to the expected GRB flux from distance r at time t. If 
we integrate this (at a time t) over all distances then we get a quan- 
tity, j r iifdr, t)dr, which is proportional to the expected GRB flux 
from that globular (at time t). The important thing about this quan- 
tity is that it takes into account the uncertainties in the reconstructed 
globular cluster and solar orbits. 

We now extend this concept to the complete set of globular 
clusters. Each cluster has a different probability per unit time of 
producing a GRB, proportional to the factor w c , defined in sec- 
tion|2.3| We can then see that the quantity 



Y(f): 



J^''=''max 



1 



-f c (r,t)dr 



(4) 



is proportional to the expected GRB flux at the Sun at time t from 
any globular cluster. In principle we integrate up to r max = oo, but 
in practice we can truncate it to a few kpc. Indeed, if there is a min- 
imum flux threshold below which the gamma ray flux is too small 
to have any significant affect on the Earth's biosphere or climate, 
then truncation is appropriate. Note that the absolute scale of ¥(/) 
is not calibrated: only relative values are meaningful. 



2.3 Weighting individual globular clusters 

Observationally, the frequency of occurrence of GRBs in individual 
globular cluster is not known. The dynamical formation of compact 
binaries, proposed progenitors of such events, is rather complex, 
involving at least two stellar encounters (see |Ivanova et al.|2008] 
|2010j >. However, the rate of GRBs in each globular cluster is ex- 
pected to be linked to the cluster properties. Several authors have 
already investigated the dependence of the compact binary forma- 
tion rate on the characteristics of the clusters. [Ivanova et al. 1 (2008 1 
found that the formation of close double neutron star binaries de- 
pends on the square of the cluster density, and that the number of 
retained neutron stars increases as the escape velocity (and thus 
cluster mass) increases. Grindlay et al. ( 2006 ) used a model where 
the formation of double neutron star binaries scales linearly with 
the neutron star number density, the velocity dispersion (and thus 
mass of the cluster) and the number of potential progenitor systems 



2 GRBs are of course discrete, rare events. Lacking information on when 
they occurred, the best we can do is to derive the probability per unit time 
of a burst for each globular cluster. 
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Figure 2. The variation of the probability density, f c (r, t), of the distance r 
between 47 Tuc and the Sun as a function of time t, shown as a grey scale. 
This scale is normalized such that the integration over r at each t is unity. 



(binaries containing one neutron star). Both models find that mas- 
sive clusters with a high concentration of stars strongly favour the 
formation of prospective GRB progenitor systems. Here we adopt 
a similar approach to these previous works and scale the expected 
GRB rate with quantities that are known for a large sample of glob- 
ular clusters. 

Specifically, assuming that GRBs are caused by neutron star 
encounters, then the GRB rate will depend on the number of neu- 
tron stars in the cluster and their encounter rate. We assume that the 
number of neutron stars scales linearly with the mass of the glob- 
ular cluster, m c , and thus linearly also with the cluster luminosity. 
The total encounter rate, T c , is given as r c oc pi (Pooley & Hut 
2006 ), where po is the central stellar number density and r, 



i.' ore 1^ the 

core radius of the globular cluster. Values for these parameters for 
our sample of clusters we obtained from Harris ( 1996 2010 edi- 
tion]^] Combining these two factors we get a quantity w c = m c T c , 
which is proportional to the frequency of gamma rays bursts in the 
clusters, and is used as the weighting factor in section [2~2| Accord- 
ingly, and as already noted in the beginning of this section, mas- 
sive clusters with high concentrations of stars at their center have a 
large GRB rate. We investigated the uncertainties of our approach 
by applying an alternative weighting scheme for individual globu- 
lar clusters. We followed Ivanova et al. ( 2008 1 and adopted weights 
proportional to p 2 , m c . With this approach we found that the typical 
uncertainties for the leading clusters is a factor of a few, with a few 
notable exceptions (see Sec. [3}. For the results in Sec.[3]we use the 
weights w c as defined earlier in this section. 

Having calculated the indivdual weights, w c , they are then 
normalised such that the sum of all weights equals 1. Here we 
used 141 clusters from Harris ( 1996, 2010 edition) where all nec- 
essary parameters are known. This, in principle, further allows us 
to estimate the expected absolute GRB rates for individual globular 
clusters by defining that a weight of 1 corresponds to the Galactic 
rate of GRBs launched in globular clusters. This galactic GRB rate 
can be calculated from the short GRB rate in the local Universe 
^Coward et al. 20 12\ and the dens ity of Milky 
l |Cole et al.|200l| >. This rate is 



of 8^ Gpc- 3 yr~' 
Way-type galaxies of 0.01 Mpc~ 
obtained for GRBs beamed towards Earth and is thus independent 
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Figure 3. As Figure[2]but for NGC 1851 



Figure 6. As Figure|2]but for M 13 
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Figure 4. As Figure[2]but for NGC 2808 



of the degree of collimation of the events. If it is assumed that the 
occurrence of short GRBs in the local Universe is dominated by 
bursts launched in globular clusters (Salva terra et al.|2008||Guetta| 
|& Stella|2009) , then the combined GRB rate of all globular clusters 
is 1(T 6 year -1 . This estimate is also consistent with the theoretically 
expected rate of short GRB production in these clusters ( |Lee et af] 
|20T0] >. 
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Figure 7. As Figure|2]but for M 15 



3 RESULTS 

Figure [8] shows the expected GRB flux, ^(f), for the case r max = 
5 kpc. This distance threshold covers 95% of all hazardous GRBs if 
a log-normal GRB luminosity distribution with log E y j so = 50.81 ± 
0.74 erg ( |Racusin et al.||20"TT| > and a critical fluence at Earth for 
a significant affect on the biosphere or climate of 10 7 erg cnT 2 
( |Melott & T homas 201 1 1 is assumed. (The profile of V P(/) has very 
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Figure 5. As Figure|2]but for Omega Cen 



Figure 8. The expected GRB flux, *P(f), at the Sun as a function of time 
before present, in arbitrary units. The vertical lines are the times of the 1 8 
mass extinction events compiled by (B ambach|2006} . 
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similar shape for other values of r max , the difference being that the 
"background" level is higher for larger values of r max , and lower 
for smaller values.) We see a significant variation. There are three 
broad peaks, at 70, 180 and 340 Myr. These correspond to times 
in the Earth's history when - within the limitations of our orbital 
reconstruction and assumptions made - we would expect a signif- 
icantly higher level of GRB flux than the average over the past 
550 Myr. 

Examining the plots of f c (r, t) for all clusters, we can identify 
those clusters which make the biggest contribution to in each 
peak: 

• Peak at 70 Myr. The main contributor is 47 Tuc, which has ten 
times the contribution to *P(7) than does the next cluster, NGC 1851 

• Peak at 180 Myr. The main contributor is again 47 Tuc, with 
several others contributing at a level 5-20 times lower, the largest 
of these being Omega Cen, M 13, and M 15. 

• Peak at 340 Myr. Once again 47 Tuc gives the largest contri- 
bution, with several others contributing at a level 7 or more times 
lower, the most significant of these being NGC 2808. 

The prominence of 47 Tuc is a consequence both of its high weight, 
w c , and the fact that it is one of the globular clusters which comes 
quite close to the Sun. All the main contributors are massive clus- 
ters that contain significant populations of dynamically formed stel- 
lar binaries. Specifically: 

• 47 Tuc has the second largest number of radio-detected msP- 
SRs (23, Ransom 20 08), dete cted by Fermi-LAT in high energy 
gamma-rays JAbdo et aL|2010) . In our weighting scheme (see Sec. 
|2,3[ > it would account for about 5% of the GRBs produced in globu- 
lar clusters. In the alternative weighting scheme (see Sec. [23} it ac- 
counts for about 1% of GRBs in globular clusters (for the following 
clusters this number is given in brackets). 47 Tuc is the dominant 
globular cluster in our study for both weighting schemes. 

• NGC 185 1 contains a msPSR in a very eccentric binary system 
with massive secondary ( |Feire et al.|2004| >. This could account for 
about 2% (1%) of GRBs from globular clusters. 

• NGC 2808 is a massive globular cluster with complex evolu- 
tionary history (Piotto et al. 2007). This could account for about 
5% (0.3%) of GRBs from globular clusters. 

• Omega Cen is the most massive globular cluster in the Galaxy, 
detected by Fermi-LAT (Abdo et al.|2010) . This could account for 
about 2% (10~ 3 %) of GRBs from globular clusters. For this glob- 
ular cluster the two different weighting schemes give the largest 
difference since it is a very massive cluster with a shallow density 
profile. 

• M 13 contains five radio-detected msPSRs (Ransom 2008). 
This could account for about 0.2% (10~ 3 %) of GRBs from glob- 
ular clusters. 

• M 15 has a double neutron-star binary that will merge within a 
Hubble time jAnderson et al.|1990) , eight radio-detected msPSRs 
( |Ransom|2008) . This could account for about 6% (2%) of GRBs 
from globular clusters. 

As mentioned earlier, GRBs are of course discrete, rare events. 
Indeed, our calculations suggest that only about 10 GRBs will have 
occurred within 5kpc of the Sun over the course of the Phanerozoic. 
Thus the true distribution of flux with time would comprise of a 
series of narrow peaks of various heights. Fig. 8 shows the expected 
flux at time (times a constant), so is the best single estimate of that 
distribution. 

By way of comparison we overplot in Figure[8]the times of 18 
mass extinction events on the Earth revealed by the fossil record, 



as compiled by (Bambach 2006). One may be tempted to draw a 
causal connection between one of these events and one of the peaks 
in ¥(?), although clearly there is a reasonable chance that one of 
these 18 events could coincide with a peak just by chance]^] It is 
nonetheless worthwhile identifying those events nearest to the three 
peaks. These are 

• Peak at 70 Myr: the famous KT extinction at 65 Myr BP, gen- 
erally accepted to have had a significant role in the demise of the 
dinosaurs; 

• Peak at 180 Myr: the late Pliensbachian/early Toarcian (early 
Jurrasic) extinction event at 179-186 Myr BP; 

• Peak at 340 Myr: the early Serpukhovian (mid Carboniferous) 
extinction event at 322-326 Myr BP, and the late Famennian (late 
Devonian) extinction event at 359-364 Myr BP. 

Whether or not a globular cluster GRB is implicated in any of these 
extinctions remains a subject for future work. 



4 OUTLOOK 

In this paper we have traced globular cluster orbits back to the be- 
ginning of the Phanerozoic eon in order to identify time intervals 
where a high flux of ionizing radiation caused by a nearby GRB 
is more likely. We found that the probability for such an event is 
far from flat with time during the Earth's history. It instead ex- 
hibits several distinct peaks, the most prominent ones being around 
70, 180 and 340 Myr BP. The main source of GRBs in all cases is 
47 Tuc. All three time intervals can in principle be associated with 
a mass extinction event, although a chance coincidence is likely. 
Therefore, to establish a link between a nearby GRB and an im- 
pact on the Earth and its biota, supporting geological signatures are 
needed. Geological signatures could comprise radiation damage of 
crystals (e.g. fossil cosmic ray tracks (Fleischer et al. 1967} or color 
shifts ( Ashbuugh 1988 1), deposition of rad ioactive isotopes {Par et| 
|al.|1998fr or elevated rates of bone cancer ( |Rothschild et al.|2003| >. 
The time intervals identified in this paper can be used as a guide- 
line to search for such signatures in the geological record. 

Finally, the current orbital parameters of globular clusters and 
the solar system are subject to considerable uncertainties. (These 
were taken into account in our analysis, and contribute to smearing 
out the probability curve.) This situation will be substantially im- 
proved in the near future with the launch of the Gaia satellite, which 
will determine the dynamics of the Galaxy with unprecedented ac- 
curacy. With better determined orbital parameters we will be able 
to constrain the past orbits more tightly, and so repeat this study to 
give results of higher confidence. 
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